function Y=lom_dis(n,Pop,U,E,I)
% Get actual Ur and LFPr from transition probas
% Get actual Ur and LFPr from transition probas
% Ordering of lambda: [jei jeu jui jue jie jiu]

%transition probas:
nei=n(:,1);
neu=n(:,2);
nui=n(:,3);
nue=n(:,4);
nie=n(:,5);
niu=n(:,6);

%hzd rates:
lambda=-log(1-n);
lei=lambda(:,1);
leu=lambda(:,2);
lui=lambda(:,3);
lue=lambda(:,4);
lie=lambda(:,5);
liu=lambda(:,6);

s=lei.*liu+lie.*leu+liu.*leu;
f=lui.*lie+liu.*lue+lie.*lue;
g=leu.*lui+lue.*lei+lui.*lei;

% I normalize Population to 1
% Pop=ones(length(g),1);

k=Pop./(s+f+g);

Uss=k.*s;
Ess=k.*f;
Iss=k.*g;

%Initialize E, U and I at their Steady-State values
U(1)=Uss(2); E(1)=Ess(2);I(1)=Iss(2);

% U=Uss;E=Ess;I=Iss;

for t=2:length(Uss)
    %get eigenvectors and eigenvalues of differential equation
    %ATTENTION: ordering is U E
    
    A=[1-nue(t)-nui(t) neu(t) niu(t); ...
        nue(t) 1-neu(t)-nei(t) nie(t);...
        nui(t) nei(t) 1-nie(t)-niu(t)];

    X(:,t)=A*[U(t-1); E(t-1); I(t-1)];
    
    U(t)=X(1,t);E(t)=X(2,t);I(t)=Pop(t)-U(t)-E(t);
end

u=U./(U+E);
l=(U+E)./(U+E+I);

Y(:,1)=u;
Y(:,2)=l;


